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Abstract 

Various flame tracking techniques are often used in hydrodynamic simulations. Their use is indispensable 
when resolving actual scale of the flame is impossible. We show that parameters defining "artificial flame" 
propagation found from model systems may yield flame velocities several times distinct from the required ones, 
due to matter expansion being ignored in the models. 

Integral effect of material expansion due to burning is incorporated into flame capturing technique (FCT) 
Khokhlov ( 1 99 5 ) | . Interpolation formula is proposed for the parameters governing flame propagation yielding 
0.2% accurate speed and width for any expansion (and at least 0.01% accurate for expansions typical in type 
la supernova explosions.) Several models with simple burning rates are studied with gas expansion included. 
Plausible performance of the technique in simulations is discussed. Its modification ensuring finite flame width 
is found. Implementation suggestions are summarized, main criterion being the scheme performance being 
insensitive to expansion parameter (thus absence of systematic errors when the burning progresses from inner 
to outer layers); in this direction promising realizations are found, leading to flame structure not changing 
while flame evolves through the whole range of densities in the white dwarf. 



1 Introduction 

In modeling astrophysical flames one is often faced with the necessity to properly track discontinuity surfaces. A 
classical situation is simulating deflagration in a near Chandrasekhar-mass white dwarf (WD) , in currently a "stan- 



dard model" of la-type supernova (SN la) phenomenon. Typical width of a flame 1 is < 1mm Timmes et al.(1992) 
which cannot be resolved, as a total size of the box for hydrodynamic simulations is of the star scale, > 10 8 cm. 

To get energetics right one needs somehow recreate flame geometry evolution without resolving responsible 
for the latter microscales. Such a problem is not untypical in physics; one hopes to mimic the processes without 
a cutoff (= crude grid cells, of size A) by cleverly prescribing velocity of the flame Df(A) and perhaps other 
parameters, making them cutoff-dependent, so that the resulting evolution coincided with (discretized) evolution 
of the exact system. At the current stage several prescriptions for D/(A) are in use. They are based on simple 
physical models of hydrodynamic instability development (essentially a Rayleigh- Taylor instability modified by 



burning; see [Khokhlov(1994) and Khokhlov(1995) , where the proposed Df(A) is also verified numerically by 



modeling turbulent burning in uniform gravity and grid resolving instability scale) , yielding Df of order of turbu 
lent velocities on scale A, St — 0.5^ Ag A, A denoting the Atwood number for fuel-ash, g being the gravitational 
acceleration (or its projection onto the flame normal). For more general recipes for Df for an arbitrary grid (giv 
ing essentially the same values for uniform Cartesian grids) see Clement(1993)| |Niemeyer fc Hillebrandt(1995) 



"renormalization invariance" of these prescriptions was numerically checked in |Reinecke et al.(2002)| , though few 
details are provided. Global performance of the numerical integration scheme at different resolutions is presented 
in detail in Gamezo et al.(2003)[ Gamezo et al.(2005)| . With Df(A) defined as a maximum of the above St (A) 



and laminar flame speed these simulations (3D with octahedral symmetry) show that the energy released at 
A m in = 5.2 x 10 5 cm (adaptive grid used) differs from that at A min = 10.4 x 10 5 cm by < 8% within the first 
2 seconds after ignition (and is practically indistinguishable before t = 1.62 s, the moment when detonation was 
triggered). Perhaps, more direct checks would be desirable (this is particularly important for showing the validity 
of adaptive-grid simulations), but it is questionable if the program outlined in the beginning of the paragraph is 
workable at all in view of intrinsic non-stationarity of the problem, instabilities and numerical noise. The cited 
works provide evidence of this invariance in statistical sense, which seems adequate in this kind of problems. All 

M.e. a region where most energy is released, and where chemical composition mostly changes from "fuel" to "ash" of a thermonuclear 
reaction of interest; usually one considers C+O burning to intermediate mass elements or, at lower densities (when quasicquilibrium 
timescale is long compared to flame passing time) to Mg, Ne, He and H. The "ash" may undergo further metamorphoses after the 
flame passes, effect of that on the flame structure and speed w.r.t. fuel neglected. 



these treatments share (in 3D, on scales A > 10 5 cm) a remarkable property of actual Df(A) being insensitive 
to microphysics, with no direct dependence on laminar flame speed. 2 

After having chosen the favorite prescription for Df[A) one starts evolving the flame surface in time; basically, 
each element of the flame is advected by the local material velocity field plus propagates with Df normal to 
itself into the fuel. One approach (level set method, LSM) is to represent a flame surface as the zero-level of a 
certain function G: flame manifold|t = {r(t) : G(r,t) = 0}; one then prescribes time evolution for G ensuring 
the needed evolution of the flame surface plus some additional physical conditions, usually taking care of G 
not developing peculiarities hindering computations. Thus at any time the flame is infinitely thin, grid cells it 
intersects contain partly fuel, partly ash. While such a picture seems physically-motivated, the exact fractions 
of fuel/ash concentrations, fluxes, etc. are computed rather crudely, the PPM-bascd hydrodynamic solvers used 
Woodward fc Colella(1984j] spread discontinuities over several grid spacings anyway. In the present paper I will 



stick to another flame-tracking prescription Khokhlov(1994)| 



This flame- capturing technique (FCT) is a "physical" way of evolving the flame 4 . On increasing dissipativc 
terms (and, actually, reintroducing them, as the only relevant in degenerate fuel burning thermal conductivity 
term is negligible any far from the discontinuity) in hydrodynamic equations one broadens the flame. Hence 
the strategy to choose such terms so that to make the flame width equal to a few grid spacings while the flame 
propagation velocity having the desired value. One can then smoothly turn the dissipative terms off far from the 
flame. In [Khokhlov( 1 995 j this idea was implemented by adding an equation for conventional "fuel" concentration 



(denoted 1 — /(r) there) to basic hydrodynamic set of equations (continuity with known sources for the species 
concentrations, momentum and energy), 

|£+u- V/ = V 2 / + $(/) (I) 

with the simplest second space-coordinate derivative diffusion-like term and step- function burning rate <&(/) ~ 
@{f — /o)(l — 0(f — I)); 5 / G [0; 1], / = corresponding to pure "fuel"; u is a local gas velocity The source 
term in the energy conservation equation then read qdf jdt, q being total heat production under burning of the 
fuel (C+O). K and $ were chosen as described above: dependence of the velocity of the flame QJ on K and 4> 
was obtained by analytically solving Q in 1 dimension, with u being constant, simply related to the Df, flame 
width was in effect estimated crudely to be of order yj K/R, its dependence on /o was not studied. This scheme 
has been used essentially verbatim by several groups since then. 

Apart from an insignificant (at the value fo = 0.3 used) error in determining Df(K, $) in Khokhlov(1995)| 
there are good reasons to study system |JQ further: 

• As described above, the technique is quite crude. The prescription for Df must vary along the flame 
front. As applied before (see thorough description in [Khokhlov(2000) ), the only factors influencing the 



Df were local gravity and mutual orientation of the front element and gravitational field. In reality the 
flame structure/speed are influenced by other factors as well, even neglecting the role of thermodynamic 
parameters variation along the front 6 : curvature of the front; independent of the curvature flame stretch, 
resulting from shear velocity gradients. Aside from these, even in ID the above procedure has flaws. Gas 
velocity u in Q will change significantly within the flame due to the gas expansion as a result of burning, 
this will make the actual flame propagation speed less than the designed Df, based on which K and $ 
were constructed as described above. We suggest that the above factors affect the Df in some universal 
way, independent of the underlying microstructure; that is model system Q must be affected in much the 



2 (Scc Khokhlov(1995) for physical discussion of self-regulation mechanism behind this phenomenon.) This, in fact, leads to low 
sensitivity of mass- burning rate to the Df(A) used as long as resolution is enough to get 2-3 generations of bubbles formed before 
substantial expansion. That was really observed in |Gamezo et al.(2003)| : multiplying St by \/2 increased released nuclear energy, 
kinetic energy and burned fraction by less than 3% for the highest resolution used (Aa; m i n = 2.6 km), but for Ai m i n = 10.4 km the 
effect was much more pronounced, about 30%, and the difference increased in time up to 1.8 sec showing that self-regulation had not 
been ac hieved. 

3 see |Von Neumann Sz Rich tmycr(1950)| for the original realization. 

4 various simulations confirm numerically basic equivalence of the two methods |Reinecke et al.(2002)| ["Niemeyer et al.(2002)| . In 
specific situations one method may have advantages over the other one. Just "by construction" no additional steps must be undertaken 
to retain sensible "flame" width in FCT; in view of the discussion in the beginning one may also hope that this "flame" response to 
hydrodynamic factors will reasonably match that of the real turbulent burning region, if implemented properly. 

5 we define 0(x) = 1 if x > 0, B{x) = if x < 0. 

6 assuming \J~A in the St definition taking care of the bulk of that dependence. 
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same way as intricately wrinkled thin flame provided the widths and velocities of resulting burning regions 
coincide on scale A. 7 

Ultimately this study would provide one with the prescription how to modify the parameters in JQ) based 
on local density, gravity, flame curvature, background velocity fields, etc. so that to match the response of 
the flame region to major hydrodynamic factors in real life. 

Apart from improving the FCT reliability this study will show the amplitude/relevance of the mentioned ef- 



fects in flame systems, a matter of constant attention [Markstc in 19641 Dursi et al.(2003) , and how different 
is the effect of these on different test systems. 

In this paper I study the effect of gas expansion on different realizations of system (QJ in ID. The higher- 
dimensionality effects require completely different methods, the results will be reported separately. In Section 
2 the problem is formulated, the numerical procedure is described. In Section 3 we qualitatively analyze the 
master equation JJJ of FCT with matter expansion for parabolic (KPP) and step-function $(/); difference in 
flame profiles helps to formulate criteria for selecting a viable realization of Q promising adequate performance 
in non-stationary simulations. Section 4 contains analytic investigation of the problem for $ being step-function, 
confirms accuracy of the numerical procedure and provides much faster way of getting the parameters; a new form 
of the diffusion term is proposed, ensuring finite flame width. Section 5 includes simple astronomically-accurate 
fits for the data obtained in Sec. 4 valid for any expansion; these are then refined to the SN la problem, we 
summarize a few "best" (in our opinion) schemes, explicit prescriptions for coding are provided. We conclude in 
Sec. 6. Appendix contains a simple analysis of performance of the scheme in real simulations. 

This work expands an earlier investigation by the author presented at FLASH Center site review (University 
of Chicago, Oct. 2004; unpublished), comments on the performance in non-steady simulations and approximate 
formulae for parameter values are added. Different from Khokhlov(1995)| form of the diffusion term (first term 
on the R.H.S. of QJ) i s a ^ so studied here for the first time. 



2 Formulation and the numerical method 

2.1 Boundary problem for finding Zfy [$(/); A] 

In this paper we restrict ourselves to a simple one dimensional model, the one where all the quantities can be 
expressed algebraically as functions of the fuel concentration field only 8 . The pressure is considered hydrostatic, 
as in deflagrations matter velocities are at least an order of magnitude smaller than the sound speed, apart maybe 
from the outer regions of the WD Timmes et al.(1992)| . The diffusion equation we consider reads 



/ being the concentration of the reaction product (1 — / that of the fuel), u local velocity, p density of the mixture 
and K the diffusivity. 

We are searching for self-similar solutions (all quantities depending on x — Dft only); as applied in FCT this 
means we neglect thermodynamical variables changing (w.r.t. flame position) on timescales of establishing such a 
steady-state burning regime. In the system of reference where the flame rests the gas velocity u{x) monotonically 

7 of course, relations like scaling oc flame width will not hold for Q-like systems as the former are critically dependent on the 
mechanism behind the structure of turbulence in gravitational field. Yet the hope is that the role of gravity is just to set scales of 
turbulent motions and hence the effective transport coefficients (like eddy viscosity/diffusivity), and knowing those alone is enough 
to model the effects concerned in the main text. 

8 this applies, say, for similar concentration and temperature fields. It is known |Zeldovich Sz Frank- Kamenetskii(1938)| that often 
used approximation when only one reaction proceeds in the front falls into this category. 

Remarks on physical nature of the constructions presented are made in several places in this paper. These are more than plain 
motivation; though the main effect - dependence of the artificial flame speed on expansion as a result of burning can be correctly 
taken into account with many bizarre variations of (all of these yielding the needed flame speed in steady-state problems, with 
further caveats like matching real u in with its models, say 131141 ) the toy systems properly modeling averaged distribution of 
density and concentration of the principal rcactant (or, usually equivalently, distribution of heat release) in real flames can fairly be 
expected to respond in the same way to the non-stationary factors ignored in this paper. We will be able to estimate only some of 
the latter in the further publications, so at any stage the naturalness of the model will remain the argument that the effects ignored 
on physical grounds are not amplified in the unphysical model beyond what physical intuition could suggest. 
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changes from — uq — Df as x — > — oo to —Df at x = +00, most of this change happening within the flame width; 
mass conservation allows to express u 

u(x)p{x) = -Dfpo, (3) 
(po being the density of the fuel far ahead of the front.) Substituting the self-similar form of solution and the p 

P 7 -Po y ' 

(enthalpy conservation; ideal gas EOS used in this section, q signifies the heat produced in the reaction) into @ 
yields a nonlinear equation for /: 

Po^ = R ff 
p dx dx 2 



^1 = K^ + Hf); (5) 



$ throughout stands for the burning rate. We mainly consider two model $: Khokhlov(1995) 



3. R for/ </<l 
1 otherwise 



and [Kolmogorov et al.(1937) 



$=iJ/(l-/) (forO</<l). (7) 
The goal is to find eigenvalues of Df for the boundary problem, JSJ) (rewritten using Q as 



[ J V 2A 

for analytic consideration) with requirements 



/(-oo) = 1 (9) 
/(+00) = (10) 
0</<l Vi€l. (11) 



For applications in the FCT A in (0 must be chosen based on the integral density jump across the flame, 



ID 



A= ^ eO^'-I). (12) 

Pfuol — Pash 

To typical WD densities pf uo i € [3 x 10 7 ;3 x 10 9 ] g/cm 3 correspond A from 0.69 to 6.3 (data adapted from 
[Khokhlov( 1995)1 , for 0.5C+0.5O composition). 

Bearing in mind the localized character of heat production/expansion under real (Arrhenius-type) burning, it 
may seem logical to improve our linear Q like 

- = 1+ A7pTT g Cf-io) (13) 
P A(l - fa) 

with p remaining unchanged in the region / < /o without burning (even though / changes with x due to 
diffusion); fuel not expanding in long diffusion tail (see below) really yields Df significantly less sensitive to A; 
this is accompanied with broader flames. Though the latter situation may seem more convenient there is no clear 
way to generalize (I13f) to more general $; moreover, as shown in the Appendix, J3J IS > f ar better motivated. Below 
we study both models (called I and II for short) to get feeling of the sensitivity of the results to the EOS. 

9 / must be continuously differentiable; left and right limits of d 2 f /dx 2 will be different wherever f = fo for the step-function 3>. 
10 A = would hold were the fuel ideal gas. Even though for generic fuel ^-(f) is nonlinear we expect the Df based on 

— = 1 + ^ with 1121 accurate enough; see Appendix for more. 
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2.2 Numerical procedure 



Upon separating <f> into a product of scale- and form-factors, $(/) = i?$o(/) onc m ay deduce the explicit 
dependence of Df (supposedly existing eigenvalues) on R and K: namely 

D f = dVKR, (14) 

where d is to be found as an eigenvalue of the boundary problem for / considered a function of a dimcnsionlcss 
coordinate 

X = x\' 

satisfying the following equation: 

and boundary conditions like ifi^l llljl . 

As l|15|) is invariant under translations in x it can generically be reduced to a first order equation by rewriting 
it in terms of 

df 
dx 

considered a function of /: 

|_,( 1 + /) + M).„. (16) 

This form is used below for qualitative and numerical analysis of the problem. Corresponding boundary conditions 
are p(0) = p{l) = 0. Eigenfunctions p{f) are nonnegative. 

To diminish the effect of singularity at p = 0, (|16ll was rewritten in terms of v — p 2 , 

I^ =< Vt?(l + //A)-#o(/). (17) 



In order not to pay extra attention to locating v — points, y/v was replaced with \/\v\ 11 . I|17|l was integrated 
by the fourth order Runge-Kutta algorithm starting from / = l,p= 0. 12 Grid spacing was constant (A/ = 1CP 5 
for most of the runs), except near the starting point and wherever a relative change in v as a result of a single 
iteration exceeded 1/16 where it was refined further. 

The d eigenvalue was then found for step-function < &o(/) by requiring p\/=i = 0. Namely, Newton-Raphson 

9(v(f=l)) 



algorithm (see, e.g. Press et al.(1992)| ) was applied, — gd — - was found simultaneously with p(f), ensuring fast 



convergence. d(A = oo) (found beforehand by solving $2'6$ ) was used as a seed at each new /o value for the first 
A in the row, for subsequent A's the previous one provided seed value for d; 4-13 iterations were enough to get d 
with 10~ 8 precision. 



3 Qualitative behavior at different burning rates 
3.1 Step-function 

We start with the model described by (JTSJ with <&o(/) = (1 — #(/ _ _ fo)- The self-similar solution 

describing flame propagation we are interested in must behave as follows: 

/ = 1 \fx < xi, df/dx (xi) = 0, < / < 1 at x > x\, lim f(x) = (18) 

x — >oo 

11 A side effect was instability of the solution next to / = 1 with respect to becoming negative with fast growing absolute value. 
To account for this v was set to whenever it became negative near / = 1. The procedure got stabilized fast, the resulting integral 
curves closely followed analytically predicted behavior, ensuring the reliability of the algorithm. The analytic asymptotes were used 
in later realizations as seed values for v near 0. 

12 It is imperative to start from / = 1 for parabolic $o(/) as a general solution for f(x) near x — > +oo (where / — > 0) is a 
superposition of two decaying exponents, and the faster decaying one is lost when integrating dp/df, thus making it impossible to 
satisfy p\ y =1 = 0. 



(moreover, it is monotonically decreasing, in accord with physical expectations; sec typical flame profiles in Fig. 1). 
Proof: as any solution of l|15(l equal to 1 at any point with non-zero df jd\ would necessarily exceed 1 nearby, 
the only other possible behavior near / = 1 would have been a solution approaching / = 1 asymptotically as 
X — * — oo. The behavior of the solution of l|15|) in that region would have coincided with that of the linearized 
equation, 

/ = ci - x/d + c-2 exp(-dx), (19) 

where 

d = d(l + l/A). (20) 

The latter does not remain in the vicinity of 1 as x ^ — oo (in fact is unbounded) for any values of the integration 
constants ci^; the only possibility for a solution to stay in [0; 1] is to become identically 1 at % smaller than some 
Xi, differentiably glued to / behaving like l|19|) as x — » Xi+- 

Similar analysis in the region where |/| <C 1 yields general solution of linearized equation 



which tends to at x 



f = ci + c 2 exp(-G?x), 
-oo iff ci = 0. 13 This completes the proof. 14 



(21) 
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Figure 1: Flame profiles (^numerically found eigenfunctions). For each /o (the ordinate of the curve intersection 
with x = 0) 3 pairs of profiles are depicted, for 1/A = 0.05, 4, and 20, larger 1/A corresponding to the curves 
intersecting x = axis at larger angles. Blue circles represent profiles for model I, red squares those for model II. 



Summing up, the desired cigenfunction asymptotically behaves like (|21|l with c\ = as x ~~ * +oo, and like 
(|T^|) as x - * Xi+ (Xi = x i s/R/K), with ci,2 such that 

/(Xi) = 1, df/d X (xi) = 0. (22) 

Presence of the three restrictions on the two-parameter set of functions (general solution of l|15|0 provides a 
hint that only under special circumstances does the required solution exist 15 . As it is shown below, this actually 
leads to a unique value for d for any fixed values of the other parameters (A and /o here), for which there exist 
solutions of the boundary problem. 

In the q = (=r- p(x) = po = const) case ()15|l actually is piecewise linear, thus the above three restrictions 
immediately yield d(A = oo,/o) as the solution of equation (cf. i|30|) : f(Xi) is then expressed in elementary 
functions as A = oo) 

f d 2 = l-e- d ' 2 . (23) 

13 Thus any solution of 1151 such that /(oo) = cannot equal zero at any finite point. It has an infinite tail, decaying exponentially. 

14 As to / monotonousncss, one can observe even more, that f(x) is convex at / < fa and concave elsewhere. Really, the 
solution of 1161 at / < /o is p = df(l + //2A), positive and increasing. If p were not monotonically decreasing in [/o; 1], one 
might consider f r , the largest / s.t. dp/df\f r = (it exists as dp/df is continuous in (/o;l) and negative near 1): p(f r ) > 0, 
d 2 p/df 2 \f r = d/A > => p(f = 1) > 0, contradicting the boundary condition. 

15 more precisely: for any fixed xi the solution f\ satisfying 1221 is unique. Another boundary condition being vanishing of / at 
X — * +oo makes the choice of \\ immaterial. That boundary condition is satisfied by one-parameter subset of solutions, and in order 
for fi to belong to that subset gencrically there must be one functional dependence among the parameters in (1151 . 



G 



One can write the solution in the limit f$ — > asymptotically as d 2 = 4-— e / fo — e J 2 f ° — ( + -\- -2- i 3 //o_ 

/o Jo Jq y ^ Jq 



To A, 



((/ O e 1//o )" 4 ); at f = 0.3 this yields d 2% smaller than the value in [Khokhlov(1995)l ; at A's of interest the dif- 
ference will be more significant. In the limit fo — > 1 d vanishes as d 2 = 2(1— /o) + |(l— /o) 2 + |(l— Jo) 3 +- • ■ (leading 
term agreeing with an estimate in |Zeldovich et al.(1985) , p. 266). 



With minor changes all said above applies to the model l|13fl with p = po V/ < fo, qualitative conclusions 
being the same. 

3.2 Burning rate (|7j) 



It is known (see Zeldovich et al.(1985)| for detailed discussion) that the spectrum of eigenvalues of the KPP 



problem is continuous, comprising all reals above some d\. In this section we show that the same holds if one 
includes the term arising from gas expansion, find similar spectrum for a wide range of 4>'s and verify these 
conclusions numerically. 

One can qualitatively analyze the spectrum for the burning rate (JJJ along the same lines as for the © described 
in the previous section. Upon linearization in the region where 1 - / < 1 (x^oo) (|15f) yields 



- - - - - d Id 2 
f=l-b+e- X+x -b-e- x ^, A± = |± VX + 1 ' (24) 

d = e?(l + l/A) as before; thus there is a one parameter subset of physical solutions, those behaving asymptotically 
like (|2U) with b+ = 0. 

By linearizing l|15(l in the region where / < 1 one gets asymptotic behavior of a solution 



d d 2 

f^b + e- x ^ + b-e- x -x, A± = -±y T -l. (25) 

There are a priori three different situations for drawing further conclusions: 

1) d < 2: any solution getting to the neighborhood of necessarily becomes negative. 

2) d > 2: any b + , 6_ > describe physically acceptable behavior, as does a certain subset of b + < 0, 6_ > 0. 
Thus one may conjecture that for all d > 2 there is an eigenfunction: it belongs to the described above 1- 
parameter subset of solutions asymptotically approaching 1 as x ~~ * ~°° an d on becoming small at larger \ it 
behaves asymptotically like 1)25(1 . exponentially approaching zero as x +oo. Essentially, this 1-parameter set 
of solutions corresponds to one of them translated arbitrarily along x> i-e- there is a unique flame profile for any 
d (the term "unique" is used below in this sense, i.e. up to translations). 

One can show that the solution cannot approach any value in (0; 1) as x — > +oo (say, by linearizing 1|15|) near 
such a value; this is of course obvious on physical grounds), and by analyzing l(TB|) one expects the / to decrease 
monotonically. Thus for the conjecture to fail, the solution of the above set ("physical" near 1) must have an 
asymptote near with unphysical b± (fo_ < in part). The numerical results below show that this does not 
happen and confirm the conjectured form of the flame profile. 16 

16 To understand better the difference between the situation here and that in the previous section one may start with some physical 
(= satisfying b.c, /(— oo) = 1) / on the left and check if it can satisfy /(+oo) = as well. As was mentioned in footnote (15) 
there is a unique solution (modulo translations) of the linearized equation with step-function <E>o, which goes to zero as x ~* 
To eliminate a possible objection that the non-zero ci makes predictions based on solutions of linearized equation questionable it is 
instructive to find a general solution going asymptotically to c\ : this reads 

/ = c1 +2a((i + ^— V^+T^-lV 1 (26) 



fo 

(up to translations in \). As ci increases from to fo J£ | increases monotonically from —fod + 1^ to 0. If for a given d the 

(unique up to translations) / going to 1 as \ ~~* ~ 00 nas ^ G [~fod(M- + 1); 0] it will go asymptotically to the corresponding ci 

x I f=fo 

as x — * namely it will be exactly 1261 where / < /o; if its derivative is more negative, it will be described by 1261 with negative 

ci until it intersectss / = at some finite x (and this is not a solution we are interested in). Easily established monotonousncss 
properties prove this way the claim that there is a unique d for which the physical solution (according to its behavior near / = 1) 
goes asymptotically to zero as x ~ * 

For the <&(/) Q, on the other hand, there are no solutions going asymptotically to any constant but as x ~ * but those 

going to represent a two-parameter set of solutions, and depending on the b± in their asymptote the ^- may take different 

x I f=fo 
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3) d = 2: (|25(l must be rewritten as f( x ) = (b + bnx) e ~" lX i and there is again a 2D domain of physical (b, bo), 
yielding / > 0, / — > as x - * +oo. In this case a meaningful burning profile exists as well. 
Analysis in this section can be straightforwardly generalized as follows: 

i) If < &o(/) goes to some (positive) constant as / — ► there are no solutions going to as x — > +oo. From the 
physical viewpoint the system reacting with finite rate in the initial state is unstable and self-similar solutions 
cannot exist. Interestingly, some problems with $o(0+) = ^ < do have eigenvalues. Corresponding cigenfunc- 
tions are identically zero on the right of some \2, f(x) ~ \ v {x ~ X2)/2) as \ — > X2~- 

ii) For <&o(/) = A*/ + °(/) at / — ► non-negative solutions going to as x — » +oo exist iff d > 2y^Z, /i > 
for positiveness of $o(/)- General solution getting to a vicinity of / = decays exponentially. Analysis near 
/ = lo/=l-/«l then suggests the following for $o = v + P>f + o(/) (below d > 2^/JI is assumed; 
/(— oo) = 1 required): 

• P > 0, pi > - for all d (the above d > 2^/Ji assumed) a unique profile exists and 3a; x : Vx < x\ f(x) = 1. 

• P > 0,/2 < - a unique profile whenever d= d(l + i) > 2^—p, (i.e. for > /i(l + 1/A) 2 the spectrum 
is additionally shrinked); again identically 1 on a half-line. 

• v = 0, p, > - a, unique solution exponentially approaching 1 as x — > — oo Vd. 

• D = fi = 0. If $o(/) = <lf r + o(f r ), r > 1, Ve? a unique solution exists. The solution of Qlfijl may be written 

as p = f f (1 + a^t - i/ r_1 + °(/ 2 + /O) . fading to / ~ ((1 - r)q X /d) ^ . If *„(/) = in some neighborhood 
of / = 1 bounded solutions exist but / — > 1 at x — * — oo from above. Not catastrophic for FCT, but unpleasant. 

• <&o ~ qf r , but r <E (0;1). / vanishes identically on the left of some xi> m its right neighborhood 
/~ ^q/2(l + r)(l-r)( X - Xl ))^. 

iii) For $o(/) = 9/ r near / = (q > 0) there are no solutions of (|16|l with p(+0) = when r < 1 - a situation 
analogous to i). When r > 1 on the other hand, there are multiple solutions going to (the ODE's peculiarity), 
hence one expects the same behavior as in ii), depending on the $o shape near / = 1. When <I>o(/) = in some 
neighborhood of / = one would expect a discrete spectrum of d, with corresponding eigenfunctions behaving 
near / = 1 according to the $o(/) near 1 as in ii). 

Summing this up, the KPP-like behavior is quite universal for $o(/) going to linearly or faster at / — ► 0; 
on the other hand, $o(f) = near / = leads to the spectrum similar to that in the previous section; positive 
$o(/) decreasing slower than linearly (or not going to zero) as / — * lead to absence of steady-state solutions. 
Profile f(x) has an infinite tail at / = 1 iff <3> 1 / ► l decreases linearly or faster. 
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Figure 2: Flame profiles at different A and d. For each c? 3 curves are depicted, for 1/A = 0.05, 5, and 20, larger 
1/A corresponding to the curves intersecting x = axis (where / was set to 0.5) at larger angles, and having 
larger widths. 

Flame profiles f(x) found numerically for $o(/) = /(l — /) are shown in Fig. 2 for four values of d. These 
seem to satisfy the boundary conditions for d > 2, whereas the profiles for d = 1 (each integrated from p(l) = 
— df /dx\f=i = 0) intersect / = at finite X with non-zero df j dx- Corresponding integral curves p(f) at the same 
values of d and A are represented in Fig. 3. Note that for d = 1 p|/=o ^ 0, whereas for d > 2 it was checked that 

values, thus making it possible to match the ^| of the solution physical near / = 1 for a range of d's. In this case /o denotes 
some intermediate value of /, between and 1. 
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by refining the grid p(0) became correspondingly closer to 0, down to p(Q) w 1CP 10 at the bulk grid spacing of 
1CP 7 (at finer grid rounding errors in real(8) dominated). This numerically confirms our qualitative conclusions. 

Note that the flame width grows fast with d. As in the original model |Kolmogorov et al.(1937)| (only asymp- 
totes near / = are essential for the argument to hold, and these have the same exponential form regardless of A) 
one may conclude that only propagation with the smallest velocity is stable. Of course, if the f(x) at some time 
corresponds to some eigenfunction above, such a profile will propagate with corresponding d. But if one considers 
a process of setting up the steady-state propagation, with initial fix) corresponding to pure fuel on a half- line 
(or / decaying with x faster than the fastest exponent A_(d = 2) = 1 in (|25Jl ) the conclusion will be that the 
resulting self-similar front will be the eigenfunction with the smallest velocity. More generally, by considering the 
evolution of some superposition ^ of the steady-state profiles found above one concludes that amplitudes of all 
of them but one will (asymptotically) decrease in time in favor of the one with the smallest d in the spectrum of 
i£ (they interact due to nonlincarity of the system) , decay rate being proportional to difference of corresponding 
A's. As a generic perturbation is a superposition of all the cigenfunctions, however small it is it will eventually 
reshuffle the profile into that for smallest d. 
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Figure 3: p(f), integrated from / = 1, p = at A and d as in Fig. 2; larger 1/A correspond to smaller p. Note 
that at d = 1 p(f = 0) is non-zero, in contrast with the d > 2 curves. 

It does not seem feasible to use KPP-like profile in FCT because of the continuous spectrum of the velocities, 
thus long times of relaxing to steady state and large widths with two exponential tails. 



4 Step-function: velocities and widths 
4.1 Analytic solution. Model I 

Let us fix the solution (thus eliminating the translational invariance) such that /(0) = fo- By going to dimen- 
sionless 

t = d X = ^fx, <f> = //A (27) 
and integrating once over £ one gets the first order equation for </>(£): 

Somewhat frivolous notation allows using this for two regions: where (f> < fa/ A (using 4>o = 1) and where 
fa/A < <f> < 1/A ($0 = 0). A is an integration constant, a priori different for these two regions. The boundary 
condition <f> = d<p/d^ = at £ = +00 fixes A = for the right half-line (£ > 0); one can further see that for model 
I left A is zero as well, by enforcing continuous differentiability of (/>(£) at £ = 0. 
One can thus evaluate the flame width as (in the original units) 

Wl . . * h + * V 1 _ M UU + *)Y 1 . (29) 



N//^l /=/0 V 2A 
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as well as the £1 (£1 = xi/d is the rightmost point where <f> = 1/A: ({1/A}) = (—00; £1] ) 

6 = -rf 2 (1 + A)- (30) 



2A, 

Note that the solution at positive £ is (in the original variables) 

/ . ai (( 1+ -).*._ 1 )- 1 .( i + ^).-*.. 

the last estimate is for |x| 3> |jjr| ln(l + 2A//o); thus for small d (achieved, say, at small A - see below) the 
region where solution decays slower than asymptotic exponential is wide (may be wider than w\ at physically 
reasonable parameter values) , and the definition of width is ambiguous (more so taking into account complicated 
f(x) dependence on the left half-line, cf. i|35l0 . Therefore we will employ another definition for comparison 

W2 = 09 (32) 

in numerical estimates. Analytically one finds (for fa > f- = 0.1) 



l + 2A//o V 2AJ 
the first term derived from the form i|31|l and the second one from the above £1. 

In the region (£1, 0) one can express general solution of i|28|) in terms of Airy (or Bessel) functions by substi- 
tution 

♦ --1 + 24, (34) 



one gets 17 



Here 



1 -(^) 1/3 (^ + ^ ln ^ + BK ^ (CT) 



(35) 



a; v*- — ) ^ 



2\ 1/2 / d'A 



and I, K are modified Bessel functions. By imposing boundary conditions 



, /o/A ata = d 2 A/6 
^ H 1/A at^l + i) 3 ^ ' (3?) 



one can get the following equation, determining eigenvalues of d: 



1 °( 1 + 1T + T-) (Wl + ^-j +K ) -I«->K = 0. (38) 



A 3a J "J \ \ 3cri 

The notation is Iq = 11/3(00), Ki = Ki/ 3 (0i), etc.; the second term differs from the first one by switching 
functions I <-+ K, i.e. it reads ^Ko(l + ...)+ Kg J (li . . . With the above definitions this equation has a unique 
solution for d for every /o, A. 



4.2 Model II 

In a similar way one may analyze the boundary problem with p given by . (D f below denotes the eigenvalue 
of (JHJ with (|9T lllf> . d - corresponding dimensionless speed, e.v. of |JSJ.) One gets exponential decay of / in (0; fo), 

17 this substitution does not introduce new integration constants, as rescaling of y (the resulting equation for is linear) leaves 
<j> unchanged. This was used in 1351 . where the coefficient at f function was set to unity. This is possible for all the solutions y except 
y = cr 1 / 3 ■ BKi/ s (a); the latter never yields eigenfunctions in our problem. 
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Figure 4: Flame velocities. The nine sequences d(A) for each model (blue circles and red squares for model I 
and II respectively) correspond (top to bottom) to fo = 0.1,. . . 0.9 with step 0.1. The worst (/o=0.1) fit 1(4*5)) for 
model I is shown. 



f = foe K x > the widths being 

Wx = K/D f , W 2 = f- f (in f- + 3 2 (l 

The equation for 

m> A ~A(l-/ ) 

in the region (£i; 0), 

dd> ~{ d>\ $ f f 

''1 + + -2i + 4^ = 0, 



i-/o 

2A 



d£ V 2/ d2 A A 
again has a general solution in terms of Airy functions 



= -1 - (tPA/4)- 1 / 3 -^- in f Ai + B Bi) (s), (39) 



ds 

which after imposing the boundary conditions as above yields the equation for d : 



(Bi^+Bi'O fAi ^-|^ +Ai j -Ai^Bi = 0. (40) 

In the above 

s = - ( 2J 2 A)- 1 / 3 ^-^(1-2/o/A) 

and subscripts at functions again represent arguments (subscript "0" means argument so = (d 2 A/4) 2 / 3 (l — 2/o/A), 
"1" similarly implies si = (J 2 A/4) 2 / 3 (l + ^Y)- 



4.3 Results for velocities and widths 

Numeric solutions of H40(l and of (|38(l arc presented in Fig. 4. 19 Figs. 5 and 6 show the wi and wi for these two 
models. Flame profiles in Fig. I were obtained by direct numerical integration of (|17fl with the d from Fig. 4. 

18 this form is more convenient here, as the argument s of Airy fun ctions Ai and Bi (b oth bounded at 0) may change sign on (£i; 0). 
In solving the equation I4UI the Airy functions normalization from |Press et al.(1992)| was used. 

19 Possibility of negative arguments for Airy functions in model II leads to infinite sequence of solutions of 141 II for the d. However, 
only the 1 st (smallest) d corresponds to a legitimate eigenfunction: the "eigenfunction" /n(x) corresponding to the n th root for d has 
n — 1 first-order poles (but it perfectly goes to as x ~ * 00 an d /(x) = IVx < Xin-) 
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Figure 5: Flame width Wi at the same A and /o- Larger widths correspond to larger /q. 



Relative difference between d's found numerically and analytically is less than 5 x 10 _3 %. It was that large for 
/o > 0.7; for /o <E [0.3; 0.6] the difference was of order 10 -8 , the accuracy required of the d in numeric procedure 
(that for solving (|38() and l|40|) was set to 4 x 10 -16 and apparently these errors played no role), and monotonically 
increased to 20 x 10~ 8 with /q further decreasing to 0.1. Quite similar numerical integration errors were observed 
for model II, reaching 8.4 x 10 _3 % (again, with bulk integration step in / being 10~ 5 ). These largest errors 
dropped to 10 _3 % when the bulk integration step was 8 times decreased. Errors in widths followed similar trends 
(and "numerical" widths were consistently larger than "analytical" ones, whereas velocities were smaller), though 
there was additional contribution from crude estimation (up to the whole grid) from integral curves x{f) obtained 
from p(f) via further trapezium integration; the discrepancy in in both models was less than 8 x 10~ 3 %. 
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Figure 6: Flame width u>2(/o, A) = — x|/ = o i j 0.9. For model I the order of /o (from larger to smaller W2 at 
1/A = 20) is 0.9, 0.1, 0.8, 0.7, 0.2, 0.6, 0.5, 0.3, 0.4 . That for model II is 0.1, 0.2, 0.3, 0.4, 0.9, 0.5, 0.8, 0.6, 0.7 . 



For the problem of immediate interest (deflagration in SN la) matter expansion is not large, thus it is worth 
trying to treat 1/A as a small parameter. Resulting asymptotic expansion looks most natural in h = l/d 2 : 



(f a -h)e^ h + h ± 



Ah + 5ti 



+ e 1/h ((fo - h f - hh 2 - 2(/ - h) (l + 



{h-hfe 2 ' h 



+ 0(A- 2 ) =0, 



from which one finds a first-order correction to the solution ho = d 2 of (|23|l . 



1 + ^| + 0(A- 2 ) ] 



-i/h 2 + hp — hpfp 
/o-e-V^o ■ 



(41) 



(42) 
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When fo is small as well this may be estimated with the aid of expansion from the end of Sec. 3.1 for ho'. 

fti = 5/o - e-V/o (J- + 1 - 6/0) - e-*/f° (Jg + 1 - 6 - 6/ ) + O(e- 3 /*>// 3 ), (43) 

this is usable up to /o ~ 0.3. Error of using does not exceed 1% for 1/A < 1 (but d = hr 1 ! 2 must be used. 
Expanding d up to 0(1/ A) leads to far worse agreement.) At 1/A = 3 the error grows to 6.7% for fo = 0.3 (1.8% 
for f = 0.2); at A > 0.69 (SN la problem) it is within 2.1% for f a = 0.3, 0.16% for f = 0.2. At these / both 
expressions for h\ give the same agreement. 

A very similar 1/A correction to (12. '{I) may be written for model II by observing that factors in front of A 
and A -1 under 1/A decomposition of l|40|l (rewritten in terms of I and K functions) coincide with the factors 
at the corresponding powers of A in the decomposition of (|38|l : therefore (|41H43|) with A switched to A provide 
corresponding estimates of d(fo,A). 



4.4 Finite width flames 



As we saw in the previous section an important drawback of KPP-like burning rates is related to exponential tails 
in the corresponding eigenfunctions f(x)- Behavior of the FCT models based on ^o(f) vanishing identically near 
/ = in non-stationary simulations is better as they admit only one steady-state solution, yet for these f(x) still 
has an infinite tail thus creating unphysically wide preheating zone in simulations hindering performance in the 
presence of large gradients. In this section we modify the diffusion term in J2J) to make the total flame width 
finite. 

Several possibilities have been investigated, 



d( P f) , d(pfu) 



dt 



dx 



d_ 

dx 



KK (f) 



'11 

dx 



(44) 



with Ko(f) = f r , r G (0; 1) seems the simplest one providing promising results. 20 With po/p = 1 + f /A this 
problem admits a unique steady-state solution with -D/ = d\/ KR, d being the eigenvalue of 



-(K Q (f) P )-d 



*o(/) 



P 







(45) 



with p(0) = p(l) = (as before p = -df/d X = -Jk/Rdf/dx). The above K (f) leads to f(x) C^smoothly 
monotonically interpolating between / = 1 Vx < Xi = — ^(1 + 1/2A) (fixing /(0) = fo; f = 1 — (x — Xi) 2 /2 + 
d(x - Xi) 3 /2 + ■•• at X -> xi+) and / - V X > X2 (f = (rd( X 2 - x)) 1/r (l + 0( X 2 - x)) at x ^ X2~)- 

Xi and the total width W3 = xi ~ Xi may be expressed in elementary functions of fo, A and d for rational 
r (or as incomplete T— function in general). Values r — 1/2 and r = 3/4 seem most adequate. Corresponding 
widths are 



in 



2 5 / 4 A 3 / 4 cT 



-4 

Arctan 



df 

--1/2 
(2/o/A) 1 / 4 
1 - (/ /2A)V2 



z , — 
— Xi = -v2A arctan 
d 



hi 



d I 

(/0/2A) 1 / 2 - (2/ /A)V 4 

(/ /2A+ 1)V2 



1 

2A 
+ d 



1+ 2A 



(46) 
(47) 



In the last expression the branch of Arctan used is Arctan € [0;7r). 

20 It might seem more physical to write the diffusion term as V (pK\7 /). This case was studied as well, normalized to unit total 
width flame profiles are more sensitive to A than the ones corresponding to 1441 : such a model is also less tractable analytically. 
The model with diffusion term of this form and with K(f) = const was also studied as a possible alternative to the original one (5) . 
Asymptotic behavior of d and the widths at small and large A is the same as presented in the previous section, hence the decision to 
stick to the original (apparently consistently performing) model. 

As with <t>o(/) we are f ree to choose from a variety of diffusion-like terms; an argument favoring the "most physical" form discussed 
in this footnote might be its more natural response to nonstationary hydrodynamic factors (following the logic of Introduction). 
However, physical transport coefficients change drastically after passing through the flame, making it questionable if putting the p 
under V is really an improvement, in view of the artificial Kq(/) dependence chosen vs effective K(f) for turbulent burning. 
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For any r in the limit A^oo 



W 3 ,r 



dr { 2Ar+l V2A/ r + 2 



O 



AY 

2 A) 



A + W^)_A(.Al_^ )+0(O 



(as in the r = case d(/o, r\A) = h 1/>2 = do(l + fx + • • ■) I/,2 > = > ^i) etc - are now somc functions of 

fo and r); as A — > W3 diverges as 



d sin irr 1 — r a Go sin 7rr 



A'- 1 (1-G 1 A) + ^ 



2 Jo" 1 
Go 1-r 



0(A). 



We used d — GqA(1 + GiA + o(A)) as in the step-function model with a standard diffusion term: d(A) dependence 
is qualitatively very similar. The fits of d(A) for r — 3/4 and various fo are presented in the next section. 
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Figure 7: Normalized to unit width flame profiles with fo = 0.6 (left) and fo = 0.2 (right). The three curves for 
each r correspond to 1/A G {0.16; 0.6; 3}, lower near / — profiles for larger 1/A. 

For a range of {fo,f} flame profiles deliver what they were expected to originally. Their width W3 may 
consistently serve as a measure of / gradients, and upon coupling (|44(l to the hydrodynamic equations one 
would really get reasonably uniform heat release within the flame width. On the contrary, for the models with 
traditional diffusion term W2 was somewhat arbitrary quantity, most of it (for larger fo and KPP to a greater 
degree) might correspond to "tails" of a profile, and the heat release in FCT would remain localized, contrary 
to the intention to more-or-less uniformly distribute it over a few cells near the modelled flame front. Somc 
representative flame profiles arc shown in Fig. 7; they are normalized to unit width (that is, reexpressed in terms 
of a new X = (x — Xi)/ (X2 — Xl)i resulting supp(<if /dX) = [0; 1]); they illustrate our "top 3 picks" to be used in 
SN la simulations: 

• r = 1/2. This has advantage that f(x) behavior near / = is the same as near / = 1, thus we are unlikely 
to introduce new problems (compared to the original r = scheme), fo = 0.2 is then convenient as the f(x) 
shape is least sensitive to A in its range of immediate interest, 1/A G [0; 2]. 

• r = 3/4, fo = 0.6. For the A's of interest corresponding flame profiles seem most symmetric overall w.r.t. 
f 1 — ► 1 — /; this is perhaps a better realization of the above idea: as the whole W3 width will be modelled by 
3-4 integration cells this approximate global symmetry seems more adequate to consider than the symmetry of 
minute regions near / = and / = 1. 

• r = 3/4, f = 0.2. Profile gradients are still uniform enough on [0.1; 0.9], they drop to zero in a regular 
manner within reasonable A% < O.25W3. At f = 0.2 the profiles seem least sensitive to A, ensuring consistent 
performance in all regions of the star. 
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5 Implementation suggestions 



Analytically found asymptotes for d(A) (Sec. 4.3) at small 1/A do not provide enough accuracy to be used for flame 
tracking in the outer layers of WD. Next order in 1/A seems marginally sufficient at /o < 0.2 yet computations 
become too involved in r ^ case (Sec. 4.4). More importantly flame capturing as presented is a general method, 
thus it is desirable to get results with larger range of validity in a ready to use form. In this section we present 
fits neatly interpolating between small and large A regions and then summarize the procedure for getting K (K) 
and R for (JJJ (respectively H44|0 in the SN la simulations. 
For model I and its modification in Sec. 4.4 

d(A) = - mi .. + n W3 /A , 9 (48) 
v ; l + m 2 /A (l + m 4 /A) 2 v ; 

with mi. ..4 obtained at each / to minimize the mean square deviation (with weights proportional to actual 
d(A)) was the simplest fit found. 21 The values shown in the graph below guarantee 0.2% accuracy for each 
fa e {0.01; 0.025(0.025)0.975} and A G [10~ 3 ; 10 5 ] studied (for any A in fact, as comparison of asymptotes 
shows). For /o > 0.3 agreement is significantly better. 
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Figure 8: Fit l|48|l parameters for r = and r = 3/4 (the last 4 curves according to the legend). Mi^ = mi^yTo. 

I have not succeeded to find a simple expression for mi...4(/o) providing good agreement for all expansions; 
this is due to delicate interplay between the two terms near 1/A = 0. 22 To overcome this a number of other 
possible fits were tested. In short, 3-parameter fits were significantly less accurate (though 3% accuracy was 
achieved, uniformly on A € (0;oo)), fits with more than 4 parameters did not yield significant improvements; 
none of accurate fits admitted simple expressions for its coefficients in terms of /o- 

This is not a major issue as /o is a parameter one can fix at some convenient value throughout the simulation. 
In |Khokhlov(1995)l f was fixed at 0.3, which essentially yielded thinnest flames throughout the A range in the 



problem 2 ^ The values for K and R used there result in actual flame speed Df = d(fo, A)v / 7o<S', reasonably close 
to S at small 1/A and /o = 0.3 (7% smaller in the WD center, however ~ 1.45 times smaller when p — 3 x 10 7 , 
0.5C+0.5O), and flame width W ~ W2(fo: A)/3Aa;; j3 — 1.5 was used, Aa; is a grid size. Now, the prescription we 
advocate for normalizes width as well as the flame speed, thus the logical criterion for /o to suggest would be the 
fastest tail (near / = 0) decay in the profile normalized to unit width wi (we discuss traditional FCT realization 

21 Different fits must be used for model II, as the asymptotic behavior of d as A — > is different. The latter was not found analytically, 

1 + -^2- + m4A~*™ 5 J 

provided 0.3% accuracy for all /o studied, and 0.45% accurate fits were achieved with the same expression when 7715 = 2 was fixed. 
fh^fh^ varied within [0.22; 0.26], this serving as an estimate for the above a. Fits in the region 1/A £ [10 5 ;10 9 ] provide strong 
evidence that V/o a = 0.25. 

22 complicated form of mi...4(/o) is related to different significance of the two asymptotic regions in the fit: for /o = 0.01 the d(A) 
becomes reasonably linear (d/A ~ const) only at 1/A > 300, whereas for /o > 0.5 this happens at 1/A > 2. 
23 A. Khokhlov, private communication. Note that this agrees with our findings, cf. Fig. 6. 
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JU here; logic is the same as in 4.4). This translates to smallest l/wzd. This, however, can be made arbitrarily 
small by taking sufficiently small / ; I would suggest to stick to fo = 0.2, as this yields \jw-id least sensitive to A 
in its interval of interest, from 0.4 to oo; normalized this way flame profiles also exhibit fairly low sensitivity to 
A. At smaller fo at small A (< 0.5) profile shapes on / E [0.1; 1] become rather nonlinear, making the choice of 
W2 as a measure for "width" more questionable. 

One can observe that (mi + mz){fo) is a significantly smoother function than both terms separately (for any 
r). This sum obviously equals da = d(/o, r |A = oo) with stated above accuracy (< 0.2%). The latter is perfectly 
well-defined (and is a solution of (|2*3|l in the r = case). d(fo, 0|A = oo) can be approximated as 

d a = [(1 - Jo 1 - 4448 exp(0.58058(l - fr 1 ) )) / /„] °' 5 , 

1 /2 

with accuracy of 0.64%; for r = 1/2 and r = 3/4 do = [2 (1 — / ) /q" 1 ] provides 1.5% accurate fits. 



fo 


r 


mi 


m 2 


m 3 


7714 


e d , 10- 4 % 


0.1 




2.9248 


0.081016 


0.23690 


0.32153 


90 


0.2 




1.9588 


0.14574 


0.26929 


0.42824 


66 


0.3 




1.4620 


0.19389 


0.32605 


0.39759 


11 


0.4 




1.1889 


0.24193 


0.30495 


0.39044 


6.1 


0.5 




1.0312 


0.29358 


0.23117 


0.40545 


11 


0.6 




0.90646 


0.34328 


0.15479 


0.42613 


7.2 


0.9 




0.44817 


0.46634 


0.015034 


0.48381 


5.3 


0.2 


3/4 


1.1683 


0.14139 


0.40205 


0.37378 


56 


0.6 


3/4 


0.81846 


0.34834 


0.14334 


0.42779 


9.0 


0.2 


1/2 


1.3969 


0.14439 


0.35064 


0.39886 


52 



Table 1: Parameters of fit Ij48|) . The last column shows the maximum relative discrepancy between the d(A) and 
its fit (with mi. .a truncated exactly as in the Table) at A g [1/3; 20], = (|Ad|/d) max . The first seven entries 
correspond to model Q), r = effectively. 

The fit parameters mi 4 for / = 0.1... 0.9, as well as for the {fo,r} values suggested in 4.4 optimized 
for A G [1/3; 50] are summarized in Table 1. The final strategy is to pick one's favorite fo (or {fo,r}) and 
corresponding mi. ..4 from the table; then at each computation step determine the local A <|12f) and needed 
Df(A) (based on local density, gravity, etc.); Ip3%)l then defines d, and (|3^|l (respectively IpT^I or l@3)) determines 
dimcnsionlcss width w (ui2 or W3). Then 

K = -£■—, R=-f/-, (49) 
d w d I w 

(where If is a desired flame width, say 3 or 4 computational cells) will lead to a flame with the desired speed Df 
and width W within the approximation adopted in this paper (see Appendix). 



6 Conclusions 

We studied a steady state ID burning governed by a diffusion equation for one species (fuel) with a source term 
representing fuel consumption; the latter depended on fuel concentration 1 — / (e [0; 1]) only. A 3D generalization 
of this system is widely used as a fiducial field governing the flame position and heat release in non-steady 3D 
simulations. Finding correct parameters of this generalization was the principal goal of the paper. 

A numerical code was written for integrating the equation and finding the flame velocity D and width W. 
The problem was reduced (by factoring out artificial diffusivity and burning rate normalization) to finding di- 
mensionless d and w depending on the burning rate shape < &o(/) & n d Atwood number A only. Several ^o's were 
studied to estimate the accuracy of the technique. The results of [Kolmogorov et al.(1937)| were generalized 
for the system with gas expansion, the spectrum of d was found continuous as in the original paper (spectrum 
dependence on the $o(/) behavior near and 1 was investigated in general as well; KPP-like velocity spectrum 
and flame profiles were found characteristic of the burning rates decreasing linearly or faster at / — > and 1); this 
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was confirmed numerically. Analytic solution for two other burning rate equation of state combinations were 
obtained, reducing the parameters finding to solving a transcendental equation (in Bessel functions). Taking bulk 
integration step of 0.001 (smaller near peculiarities) in direct integration (in /) was enough to ensure coincidence 
of the results with these analytic ones within 1% for both d and w. 

Asymptotic behavior of d and w at A near 1 (no expansion) and co (large expansion) was established. 
Interpolating formulae were proposed for these quantities valid for any expansion coefficient, thus providing an 
efficient method for finding the parameters for the equation describing flame propagation in 3D simulations. 
Another variation of the scheme was proposed and studied, which reslulted in finite flame width, hence promising 
better performance in non-stationary simulations. Specific realizations were pointed out, effectively uniformizing 
flame profiles (further improvement beyond widths normalization, advocated for with any variation of Q), thus 
ensuring the same performance at all stages of simulation. 

Brief analysis showed that the flame shall propagate with the prescribed velocity Df and shall have the 
prescribed width W with high accuracy, being insensitive to details of the actual burning rate as long as most of 
the energy is released within its width and as long as the matter enthalpy scaling H oc 1/p holds on flame width 
scale. This is the case in the SN la problem, one usually defines a "flame" for the first condition to hold, and the 
second one is true for the first stage of deflagration in a white dwarf (while all the effectively coupled components 
of the matter can be considered as ultrarelativistic gas), in which virtually all of the nuclear energy is released. 

We argued that the restriction of burning rate depending on / only is a good approximation; burning inside 
the flame in most simulations is described by accurately modeling burning out of one element, while effect of 
burning of other species estimated crudely. This is physical — if several reactions have comparable rates, they 
can be effectively described by one / variable; if one reaction is significantly faster, the flame is influenced by it 
to a higher degree, and the approximate estimate is enough. 

For the boundary case, when different reactions may be important in different flame zones (or when Hp 
changes significantly within the flame) actual velocity of the artificial flame speed may differ from the value 
based on which the FCT was modeled. One might then try to consider steady state solutions of case-specific 
generalizations of with several species and temperature involved. As another class of effects caused by finite 
flame width are effects of flame curvature and matter velocity gradients on the flame propagation. These will be 
a subject of future publications, as will be testing of the results. 



Appendix 

Here I present some expectations of the FCT performance with the new parameters in real- world simulations. 
Basic setup is the following: the physical system consists of some reactant, its concentration denoted g(r, t), u, p 
denote its velocity and density fields. The burning is modeled as described in Introduction. First consider a 
steady state burning; in such the quantities are distributed according to {JJ an d (ID without loss of generality) 

dg_ = G(H 7 g) 
dx PoDf 

— =q^ (51) 
dx dx ' 

continuity equation and equation of state having been used to eliminate u{x) and express burning rate of g as 
a function of enthalpy H instead of the temperature (considering the pressure constant H and g are the only 
thcrmodynamically independent quantities. Equation for g effectively decouples because of Ij51|l form, dictated 
by FCT. g is generally a vector, describing concentrations of all the species modelled. In some realizations 
instead of l|50|l one simply determines g as some pre-defined function of /, that making the analysis here more 
straightforward) . 

A crucial observation at this point is that H = co ™ st would be enough to get (from JUJ) an expression for u 
in (J2J coincident with that used in model I of Sec. 2, thus tautologically yielding the same D (and then the same 
W). Rest of the equations would not matter: given / l|51|l would determine H(x), while l|50[) would define g(x). 
The boundary conditions for H and g are satisfied automatically: gradient of H tends to zero at infinities; that of 
g behaves the same at +00, assuming as usual fast reaction rate decay at low temperatures; and g{— 00) = for 
usually considered G(H,g) cx g m ,m > 2 and G(H (— co), g(+oc)) not very small (otherwise no reaction proceeds). 
On the other hand, different H(p) would yield different u(x), and the resulting eigenvalue for D might differ 
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drastically. As a ready example consider model II of Sec. 2 for prescribing K, R. Suppose in reality H does 
scale oc l/p. One can check that the actual flame speed will be d\i/d\ times smaller (that is 7.6 times smaller 
for /o = 0.9, 1/A = 20); one can expect comparable discrepancy for model I used for modeling a system with 
d^/dH\ p=const being much smaller at H(+oo) (cold fuel) than at H(—oo) (hot ash). Thus, the distribution of 
density within the flame is important, not just an integral jump. 

Fortunately for the method, H = for ideal (=non-interacting in statistical physics sense, degeneracy and 
statistics do not matter) ultrarelativistic gas, and the matter in WD core can reasonably be considered as such. 
It remains ultrarelativistic during the burning, thus making assumptions behind model I applicability fulfilled. 24 
There are no clear preferences as to what value for /o to use at this stage. 

In non-stationary burning, © (which would have yielded the prescribed flame velocity, time-dependent in this 
case) cannot be derived from (Q, u will explicitly depend on coordinates and time. Far from caustics though one 
may hope that dependence to be a small correction to |u(/)| from the previous paragraph. Changing cross-section 
of streamlines (due to flame curvature and stretch), pressure changes in the gravitational field (and thus violation 
of Hp = const relation 25 ) are examples of this phenomenon; quasi-steady state burning is possible in this case, 
and the flame speed of such is a likely real flame speed in the coupled system of hydrodynamic equations and Q . 
Linearizing the perturbed u near the flame position one might find the effect of nonstationarity being possible 
to model with some equivalent flame with curvature. We found the geometric effects not very pronounced, being 
< 5% for flames with radius of curvature 10 times their width and above. However this whole discussion is super- 
ficial in that conditions far from the flame may still have a strong effect on its speed, as the boundary conditions 
are at infinities (however the models with finite width, like those studied in Sec. 4.4, must be more robust). This 
was seen in Section 3, similar complications occur in the problems with flame curvature. Real field-testing of the 
flame evolution is needed for the final verdict. 
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